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ABSTRACT 

Feedback from super-massive black holes (SMBHs) is thought to play a key role in reg¬ 
ulating the growth of host galaxies. Cosmological and galaxy formation simulations 
using smoothed particle hydrodynamics (SPH), which usually use a fixed mass for 
SPH particles, often employ the same sub-grid Active galactic nuclei (AGN) feedback 
prescription across a range of resolutions. It is thus important to ask how the impact 
of the simulated AGN feedback on a galaxy changes when only the numerical reso¬ 
lution (the SPH particle mass) changes. We present a suite of simulations modelling 
the interaction of an AGN outflow with the ambient turbulent and clumpy interstellar 
medium (ISM) in the inner part of the host galaxy at a range of mass resolutions. We 
find that, with other things being equal, degrading the resolution leads to feedback 
becoming more efficient at clearing out all gas in its path. For the simulations pre¬ 
sented here, the difference in the mass of the gas ejected by AGN feedback varies by 
more than a factor of ten between our highest and lowest resolution simulations. This 
happens because feedback-resistant high density clumps are washed out at low effec¬ 
tive resolutions. We also find that changes in numerical resolution lead to undesirable 
artefacts in how the AGN feedback affects the AGN immediate environment. 
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1 INTRODUCTION 

Feedback from AGN is often invoke d in galaxy forma¬ 


tion and cosmologi c al simulations (e.g., ISoringel et a U 


Schave et al. 2010; IPubois et al. 1 1201 ISchave et all 


2005 


2015 


Vogelsbe rger et al.ll2014l) as well as in semi-analytical models 
(e.g., [Bower et al .1120061 : ICroton et al.l 1200(1 iFanidakis et al.l 
l2012h in order to quench star formation in galaxies at 
the high-mass end of the mass function and reproduce a 
number o f observational correlation s such as the Mbh — cr 
relation (iFerrarese fe Me rrittl l20 00l: iGebhardt et"aI1 l2000l : 


iTremaine et al.l 120021 : iKormendv &; Holl2013h . The general 
premise in such models is that the AGN provide a source 
of negative feedback, clearing gas from the host galaxy and 
inhibiting further star formation and AGN activity. 

Outflo ws on kpc scales with velocities ^ 1000 km 


(e.g., Cano-Dfaz et a bl 120121: iMaiolino et al. I I 20 I 2 : 


ICicone et, al. 1 12014 I2015I: iTombesi et al JI2015D and momen¬ 
tum fluxes exceeding the radiative output of the AGN , 
Pagn — Tagn/c, by factors of upto ^ 30 (IPunn et al. 120101: 
Feruglio et al. I2OI0I: iBautista et al. 120101: lllupke &; Veilleuxl 


20 Ilf Is tumm^tjd. I 12011I: iFaucher-G-iguere et al 


PmcherTliguere V Quataert] 20121 : iGenzel et al 


2012 


20141 


to be driven by AGN. Such observations provide compelling 
evidence that AGN can indeed have an impact on the 
host galaxy, playing an important role in establishing 
observed correlations and thus vindicating the use of AGN 
feedback in simulations and semi-analytic models (see also 
McNamara &; Nulsenl 120071 : iFabianl 120121 : king &; Poundsl 
201,4) . 

Observations of local (z £ 0.1) AGN have found that 
~ 40% of systems ho st “ultra-fast” outflow s (UFOs), with 
velocities of v ~ 0.1 c ([Tombesi et al.ll2010al H) at small radii. 
Typically such outflows have mass outflow rates M ou t ~ 
0.1 M© yr -1 a nd kinetic en ergy fluxes M out v 2 /2 0.05LEdd- 
Models rtKingl 11200 3, l2005h show that when these outflows 
impact upon the ISM, the wind shock can reach tempera¬ 
tures of order ~ 10 10 — 10 11 K. When radiative cooling of 
the wind is inefficient, it expands adiabatically and has the 
potential to drive the high velocity outflows discussed above 
and clear out significant fractions of gas from the host galax y 


and clear out significant tractions ot gas trom the host galaxy 
(|Faucher-Gigubre &; Quataertll2012l : IZubovas &; King! 20121 ) . 


Pespite the success of cosmological simulations in re¬ 
producing large scale observations (e . g.. ISchave et al.l 2010j; 
McCarth y et al.l 120101: Fabian et al.l l2010l; IPlanelles et al l 


Tombesi et al.ll2015l ) have~been observed and are believed 2013 : Vogels berger et al.ll2014l : ISchave et al.ll2015h . they are 
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unable to resolve scales small enough to probe the “AGN- 
engine” and thus provide limited i nsight into th e exact pro¬ 
cesses driving AG N feedback, see ISchave et al ] (1201 Fill and 
ICrain et, al.1 ((201 5) for a detailed discussion. Therefore sim- 
ulations only model the effects of the feedback on the ISM, 
as opposed to the feedback mechanism itself. Typically such 
models have to be tuned , that is, free parameters of the 
feedback and other prescriptions have to be varied until a 
reasonable fit to a set of calibrating observations is found. 

This unfortunate situation is unlikely to be drastically 
improved any time soon because the numerical and physical 
modelling challenges in AGN and star formation feedback in 
cosmological simulations are so great. Nevertheless, in the 
interests of the field, it is only fair to ask the question: does 
this approach create numerical artefacts that may influence 
predictions of the simulations in a systematic way ? 

To give an example, consider how the SMBH mass, 
Mbh, can be limited by a feedback argument. Suppose 
that our model for SMBH feedback contains a parameter 
cbh = £f€r that defines the fraction of SMBH rest mass en¬ 
ergy that goes into the AGN outflow, cbhMbhc 2 , where e r 
is the radiative efficiency of the black hole and ef is the effi¬ 
ciency with which the radiation couples to the surrounding 
gas. Some of this energy may be lost in the outflow-ISM in¬ 
teraction, for example to radiation in cooling shocks or by 
escaping the galaxy through low density voids (see below), 
so effectively only a fraction, cism, of the feedback energy 
impacts the host galaxy gas. In this scenario, the maximum 
SMBH mass is then limited by 

CISM€BhMbhC 2 = MgasCT 2 , (1) 


where M gas is the mass of the gas in the host galaxy that 
AGN feedback needs to remove from the galaxy and a is 
the ID velocity dispersion. From this simple analytical ar¬ 
gument the black holes mass should be determined by the 
efficiency parameters suc h that Mbh oc (^ism^bh) -1 - A sim¬ 
ilar conclusion is found bv lBooth fc Schavei (12010 ) who show 
that Mbh oc (cbh) -1 , where cbh is a free parameter of their 
feedback model. 

AGN feedback is often implemented in galaxy forma¬ 
tion simulations as a sub-grid model for which the black hole 
efficiency parameter, cbh, is set by hand. cbh is often cali¬ 
brated in order to rep roduced the observed l ocal black hole 
scali ng relations fe.g. iDi Matteo et al.l 20051: Soringel et al.l 

^21 


l200fil:ISiiack7 et al 1l2007l:lBoot,h~ & Schavell2009l ) . with typical 

values of e r = 0.1 and ef = 0.05 — 0.15. However, cism, which 
cannot be directly set by the simulator, is governed by the 
ISM modelling i.e. details of the hydrodynamics, radiative 
cooling and any other sub-grid ISM routines used in the sim¬ 
ulation. This provides an explanation as to why values for ef 
can di ffer between simulations. As noted in lBooth fe Schavei 
(l2009h their value of ef — 0.15 d i ffers from the value of 
€f = 0.05 used by ISpringel et al.l ( 2005 ) due to compen¬ 
sating for differences between sub-grid IS M modelling. This 
sugges ts that the eff ective eisM is smaller in lBooth fc Schavei 
(|2QQ9h compared to lSoringel et al.l ([20051 ) . 

From the arguments above, when a simulation is com¬ 
pared to observations, a constraint is obtained not on cbh 
directly but on the product cismcbh• The danger here is 
that cism is dependent on the numerics and hence the value 
obtained for cbh when calibrating simulations against ob¬ 
served black hole scaling relations does not actually di¬ 


rec tly tell us about th e AGN physics (as already discussed 
by ISchave et al.l 12015l ). We note that e/ is never consid¬ 
ered a prediction of the subgrid AGN feedback models and 
that the only requirement for self-regulation of the SMBH 
growth to occur is that e/ is non-z ero. Further , it is i nter- 
esting to note that both the OWLS ([Schave et al.ll201Qh and 
EAGLE (Sc have et alJ 120151 ) cosmological simulations had 
large differences in resolution and subgrid physics, but used 
the same value of ef = 0.15. This choice did however re¬ 
quire an increase in the temperature increment of particles 
heated by AGN feedback in higher resolution simulations 
([Schave et al.ll2015l : ICrain et al.ll2015h . This parameter is set 
by hand and effectively controls the value of cism- The inti¬ 
mate relationship between e/ and cism, evidenced by these 
large-scale simulations, shows that it is important to un¬ 
derstand any potential numerical trends in cism, for exam¬ 
ple with resolution, before drawing conclusions about AGN 
feedback mechanisms. Investigation of these mechanisms is 
a logical next step in galaxy evolution simulations. 

In this paper we perform a resolution study in order 
to better understand how numerical resolution can affect 
the cou pli ng between t he SMBH feedback and the ISM. As 
in dBourne et al.l [20141 . hereafter BNH14), to achieve a cer¬ 
tain degree of realism in modelling the clumpy ISM of real 
galaxies, we impose a turbulent velocity field upon the ini¬ 
tial smooth gas distribution, and allow clumpy structures 
to develop before they are hit with the SMBH outflow. We 
vary SPH mass resolution over four orders of magnitude, 
and we also vary the SMBH feedback implementation and 
the cooling prescription used in order to minimise numerical 
artefacts. Our numerical simulations allow us to test whether 
there are numerical trends in cism for a single SMBH feed¬ 
back event. Briefly, our main conclusion is that, in the sce¬ 
nario studied, feedback in low resolution simulations is far 
more effective at destroying galaxies than it is in higher reso¬ 
lution simulations. This indicates, at least qualitatively, that 
cism is resolution dependent. 

The paper is structured as follows; section 2 outlines 
the numerical method and how the simulations are set up, 
section 3 highlights the results of the simulations, section 
4 discusses the implications of these results, both physical 
and computational, and finally in section 5 we summarise 
the outcome of this work. 


2 SIMULATION SET-UP 
2.1 Numerical method 

We implem ent the SPHS0 formalism a s desc ribed in 
iRead et ahl (|2010h and iRead fe Havfieldl (|2012h . within 
a modified version of the N-body/hydrodynamical code 
GAD GET-3 , an updated version of the code presented 
i n ISpringej (120051). The second-ord er Wendland kernel 
(|Wendlandlll995l : iDehnen fe Alvlliorj ) is employed for both 
SPH calculations (using 100 neighbours), and weighting of 
the AGN feedback. The simulations are run in a static 
isothermal potential with a mass profile which follows 


1 Smooth Particle Hydrodynamics with a high-order dissipation 
Switch. 
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M(R) = —R = 

CL Gr 


( 2 ) 


where M a = 9.35 x 10 9 M 0 is the mass within a radius 
of a = 1 kpc and cr po t = \JGM a /2a ~ 142 km s _1 is the one 
dimensional velocity dispersion of the potential. In order to 
prevent gravitational forces diverging at small radii we apply 
a softening length of 0.1 pc. 

An ideal gas is used for all simulations with gas pres¬ 
sure given by P = (7 — 1 )pu, where p and u are the den¬ 
sity and internal energy per unit mass of the gas respec¬ 
tively and 7 = 5/3 is the adiabatic index. Note that the 
mean molecular weight for the gas is calculated self consis¬ 
tently in our simulations, however for simplicity we assume 
p = 0.63 when plotting temperature. In our fiducial runs, 
for gas temperatures above T = 10 4 K, we use a modified 
version of the optic ally thin radiative cooling function of 
ISazonov et al.1 ( 2005 b which includes Bremsstrahlung losses, 
photoionisation heating, line and recombination continuum 
cooling and Compton heating and cooling in the presence 
of an AGN radiation field. For comparison we also carry 
out runs using the same prescription but neglect the effect 
of inverse Compton (IC) cooling against the AGN radia- 
tion field. This is in light of recent theoretical predictions 
(|Faucher ^fflgu|£e_^_QMtaeij I 2 Q 12 ) and observational con¬ 
straints ( Bourne fe Navakshinll2013 i that suggest UFOs are 
always energy conserving a nd do not c ool via IC processes 
as was previously believed llKind 120038). Be low T = 10 4 K, 
cooling is modelled as in lMashchenko et al.l (120081 1. proceed- 
ing through fine structure and metastable lines of C, N, O, 
Fe, S and Si. For simplicity, solar metalicity is assumed for 
all cooling functions. 

We impose a ‘dynamic’ temperature floor such that gas 
cannot cool below a temperature of 


Tfloor = P 1 ^ 3 ~ —77— (AUgb^SPH) 2 / 3 

7T/CB 

- 350 ( p - ] 1/3 (-*-) (-EH*-') ^X 
“ \ 10 - 22 gcm -3 ) VO.63/ Vl 6 OOM 0 ; 

(3) 

where p and mspH are the density and mass of an SPH 
particle, respectively, and iV ng b = 100 is the typical num¬ 
ber of neighbours. Such a temperature floor manifests itself 
as a polytropic equation of state with an effective 7 = 4/3 
and is used for purely numerical reasons to guarantee that 
the Jeans mass is independent of density and Jeans length 
scales with the SPH kernel smoothing length. This ensures 
that gas clouds are able to collaps e while avoiding spuri¬ 
ous f r agmentation due to resolution ([Robertson fe Kravtsovl 
120081 : ISchave fc Dalla Vecchiall2008l b This method, or vari¬ 
ants upon it are widely used in galaxy formation and cosmo¬ 
logical simulations alike (e.g.. ISchave &; Dalla Vecchial [2008: 
iHobbs et al.l [2013) and thus, despite not being physically 
motivated, is an important ingredient in our study if we are 
to compare to resolutions similar to those achieved in cos¬ 
mological simulations. 

SPH particles that have reached the temperature floor 
and have a density above p — 10 -22 g cm -3 are considered 
star forming. The properties of the temperature floor en¬ 
sures star formation follows a Jeans instability criterion. We 
employ a probabilistic approach to convert a fraction of this 


gas into stars. Similar in fashion to iKatzl (|l992h . the proba¬ 
bility of a SPH particle being converted into a star particle 
in a given time step At is given by 

P = 1 - exp ( 4 ) 

where csf = 0.1 is the assumed star formation efficiency and 
Tff ~ ^3ir/32Gp is the local free-fall time of the gas. Newly 
formed star particles are equal in mass to the SPH particles 
and only interact with other particles through gravity. 


2.2 Initial conditions 

The simulations presented here follow a similar setup to 
those presented in BNH14. We wish to investigate the im¬ 
pact of AGN outflows on ambient gas in the host galaxy un¬ 
der realistic conditions, which should certainly include the 
fact that the ISM is very non-homogeneous, that is, clumpy. 
To achieve that condition in the controlled environ ment of 
an isolated simulation, similar to lHobbs et al.l (|201ll b a tur¬ 
bulent velocity field is imposed upon the gas. A sphere of gas 
(cut from a relaxed, glass-like configuration) is seeded with a 
turbu lent velocity fiel d using the method of lDubinski et al.l 
(Il995l h as described in I Hobbs et al.l ([201 lh . Assuming a Kol¬ 
mogorov power spectrum with P v (k) ~ /c -11 ^ 3 , where k is 
the wavenumber, the gas velocity can be defined as v = 
V x A where A is a vector potential whose power spectrum 
is described by a power-law with a cutoff at k m in — l/-R ou t 
(where R out is the outer radius of the system and defines 
the largest scale, A max = 27r//c m i n , on which turbulence is 
driven). The velocity field is generated by sampling A in 
Fourier space. At each point (k x ,k y ,k z ) the amplitudes of 
the components of Ak are drawn from a Rayleigh distribu¬ 
tion with a variance given by < |Ak | 2 > and phase angles 
distributed uniformly between 0 and 2i r are assigned. The 
last step is to take the Fourier transform of v^ — ikx Ak in 
order to obtain the velocity field in real space. 

The desired set up for the gas distribution, on which 
the AGN feedback acts upon, consists of a clumpy gaseous 
shell with a radial range from 0.1 kpc to 1 kpc and a 
10 8 Mq black hole at the center. This is achieved by first 
setting up a gas distribution which initially follows a sin¬ 
gular isothermal sphere potential from — 0.1 kpc to 
J^out 8 1 kpc. The gas mass fraction within this shell is 
/ g = Mg/Mtotai = 0.16, giving a total initial gas mass 
M g 1.6 x 10 9 M©. The system, which initially only consists 
of SPH particles and the central sink particle, is then allowed 
to evolve under the action of a turbulent velocity field for 
1 Myr, resulting in a clumpy gas distribution. The turbu¬ 
lent velocity is normalised such that the root-mean-square 
velocity, fturb — cr — 154 km s -1 and the gas temperature 
is initially set to T ~ 5.6 x 10 5 K, such that the shell is 
virialised. 

The black hole is modelled as a 10 8 M 0 sink particle. 
During the relaxation period gas is added to the sink particle 
if it falls within our desired inner boundary for the initial 
condition of 100 pc. At the end of the relaxation period the 
sink particle mass is reset to our desired black hole mass of 
1O 8 M 0 and the accretion radius is set to 10 pc. This results 
in particles at small radii with prohibitively small time steps 
being removed whilst allowing us to still be able to follow 
the inflow of dense filaments to small radii during and after 
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the AGN outburst. However to prevent the removal of gas 
directly heated by the AGN feedback, SPH particles that 
are not bound to the collective mass of the sink particle 
and background potential (within the SPH particles radial 
position) are not accreted. Here we present simulations that 
initially have iVsPH = 10 3 ,10 4 ,10 5 and 10 6 SPH particles. 


2.3 AGN feedback model 


Even at the resolutions presented in this paper we are 
unable to directly model the feedback mechanism of the 
AGN, however we can model the effect of the feedback 
on the ISM. Models of UFOs colliding with the ISM have 
been particularl y successful i n exp l aining observational cor¬ 
relations (e.g.. | Kind 20031. 120051: IZubovas <U Kind 120121 : 
iFaucher-Giguere fe Quataertl l2012 h In these models the 
UFO, with a velocity v ~ 0.1c, shocks against the ISM, driv¬ 
ing a reverse wind shock and a forward shock in the ISM. 
The wind shock can reach temperatures of ~ 10 10 — 10 11 K 
and expand through thermal pressure, driving out material 


of the ISM. As in iCosta et al 


it is the effect of the 


reverse wind shock that we attempt to m i mic in our feed¬ 
back method. Similar to iDi Matteo et al.l (|2005h we ther¬ 
mally couple the feedback to neighbouring gas particles in 
a kernel-weighted fashion. During a time step of length At, 
the energy released by the AGN is given by 


F/therm — CfF/AGNAt 


(5) 


where ef = 0.05 is the efficiency with which the AGN lu¬ 
minosity couples to the ambient gas, as defined in the in¬ 
troduction and Lagn is the AGN luminosity. Our chosen 
value for ef is physically motivated by models of UFOs 
which are expected to have a ki netic luminosity Ekin.UFO = 
(e r /2 )LAGN — 0.05 Lagn (e.g., lKin3l2005l : [Zubovas & Kind 
l2012l l . For simplicity we set the AGN duration to 1 Myr and 
Fagn to the Eddington luminosity, 


FEdd — 


AttGMbuc 


(6) 


where Mbh = 10 8 M© is the black hole mass and k = <tt /m p 
is the electron scattering opacity (where <tt is the Thompson 
cross-section and m p is the proton rest mass) and G is the 
gravitational constant. The energy given to an SPH particle, 
Fdnj, is then given by 


F/inj 5 k — F/therm 


ttisphIU( rk — ^bh, hbh) 

Pg(r B H) 


(7) 


where mspH is the mass of an SPH particle, W(r^ — 
rsH, Fbh) is the kernel weight of the SPH particle relative to 
the black hole, Fbh is the black hole smoothing length, cal¬ 
culated over fBH.N n gb neighbours (see table [I]) and p g (r*bh) 
is the gas density at the location of the black hole. This ap¬ 
proach ensures that gas closer to the black hole is heated to 
a higher temperature than gas further away. The total mass 
heated per time step is given by 


A/heat — /BHiVngbmsPH , 


( 8 ) 


where /bh is the ratio of the number of black hole neigh¬ 
bours be heated. We consider two main scenarios, one in 
which /bh = 1 at all resolutions and one in which we ap¬ 
proximately heat a fixed mass at each resolution and so set 
/bh — Asph/10 4 . This choice of /bh is a balance between 



Figure 1. Gas density distribution at 1 Myr for simulations with 
10 3 (red), 10 4 (black), 10 5 (green) and 10 6 (blue) particles. Both 
the high and the low density tails of the gas density distribution 
are better resolved as the mass resolution of the simulation is 
improved. 


heating a sufficient number of particles in the lowest reso¬ 
lution simulations and not heating an excessive number of 
particles at high resolution. 


2.4 Summary of simulations 

A summary of the simulations is given in table [T] We use a 
nomenclature of the form FNYZ or FMYZ, where “FN” sig¬ 
nifies that a fixed number of SMBH SPH particle neighbours 
are heated by the feedback independently of the SPH parti¬ 
cle number used in the simulation. This means that /bh = 1 
for such simulations. “FM”, on the other hand, stands for 
a fixed mass of SMBH neighbour particles being heated. 
In these runs the number of SPH particle neighbours over 
which the SMBH feedback is spread depends on the numeri¬ 
cal resolution of the simulation, and we set /bh = Ysph/10 4 , 
at all resolutions. The number Y= log 10 (iVspH) encodes the 
total number of SPH particles used. Finally, Z is either “h” 
or “c”, and marks runs in which cooling due to IC processes 
is included (c) or not (h). 


3 RESULTS 

3.1 Pre-feedback properties of the ISM 

Before investigating how the feedback interacts with the ISM 
we compare the properties of the ISM itself at different res¬ 
olutions. Figure [T] shows the density distribution for the gas, 
at different resolutions, after 1 Myr i.e. just before the feed¬ 
back turns on. At this point in the simulation the gas dis¬ 
tribution is identical for all of the runs at the same reso¬ 
lution, e.g. the blue curve in the figure is the same for the 
runs FN6c, FN6h, FM6c and FM6h. Figure [T| shows that 
the lowest resolution runs, with 10 3 particles, probe a much 


2 Compton heating due to the AGN radiation field is included 
for gas with T ^ 1.9 X 10 7 K in all simulations. 







































The resolution bias 5 


Run 

Nsph 

m SPH (M©) 

/BHN ngb 

cooling 


FN3c 

10 3 

1.6 X 

10® 

10 2 

Sazonov et al. 

(20051 

FN4c 

10 4 

1.6 X 

10® 

10 2 

Sazonov et al. 

(20051 

FN5c 

10® 

1.6 X 

10 4 

10 2 

Sazonov et al. 

(20051 

FN6c 

10 6 

1.6 X 

10 3 

10 2 

Sazonov et al. 

(20051 

FN3h 

10 3 

1.6 x 

10® 

10 2 

Sazonov et al. (2005b no 

Compton cooling 

FN4h 

10 4 

1.6 x 

10® 

10 2 

Sazonov et al. (2005b no 

Compton cooling 

FN5h 

10® 

1.6 x 

10 4 

10 2 

Sazonov et al. (2005b no omDton cooline 

FN6h 

10 6 

1.6 x 

10 3 

10 2 

Sazonov et al. (2005), no 

Compton cooling 

FM3c 

10 3 

1.6 x 

10® 

10 

Sazonov et al. 

(2005) 

FM4c 

10 4 

1.6 x 

10® 

10 2 

Sazonov et al. 

(2005) 

FM5c 

10® 

1.6 x 

10 4 

10 3 

Sazonov et al. 

(2005) 

FM6c 

10 6 

1.6 x 

10 3 

10 4 

Sazonov et al. 

(2005) 

FM3h 

10 3 

1.6 x 

10® 

10 

Sazonov et al. (2005b no 

Compton cooling 

FM4h 

10 4 

1.6 x 

10® 

10 2 

Sazonov et al. (2005b no 

Compton cooling 

FM5h 

10® 

1.6 x 

10 4 

10 3 

Sazonov et al. (2005b no 

Compton cooling 

FM6h 

10® 

1.6 x 

10 3 

10 4 

Sazonov et al. (2005b no 

Compton cooling 


Table 1. Summary of simulations showing (1-r) run name, initial number of SPH particles (A/sph), mass of a single SPH particle (?71sph), 
number of black hole neighbours heated during feedback (/BH-^ngb) an d the cooling prescription used. Run nomenclature takes the form 
F XYZ where X defines whether the thermal energy of the AGN feedback is deposited into a fixed number of neighbours ( N) or a fixed 
mass (M) at all resolutions, Y— log 10 (iVspH) and Z defines runs in which cooling due to IC processes is (c) and is not (/i) included. 



Figure 2. Gas temperature distribution at 1 Myr for simulations 
with 10 3 (red), 10 4 (black), 10 5 (green) and 10 6 (blue) particles. 
The low-temperature (dense) gas is completely unresolved in the 
low resolution simulations. 

narrower density range than the runs in which 10 6 parti¬ 
cles are used. The highest resolution runs thus resolve the 
density distribution tails at both the low and high density 
ends. This means that with improved resolution we are able 
to better distinguish the high and low density phases of the 
ISM, which, as we show below, can have a large impact on 
the efficiency of AGN feedback. 

Figure [2] shows the temperature distribution for the gas, 
at different resolutions, after 1 Myr. In this figure we can see 
that in the lowest resolution simulations (10 3 particles) the 
gas remains warm, T > 10 4 K. This is due to the dynami¬ 
cal temperature floor that we employ (which is a common 
approach in cosmological and galaxy formation simulations, 
see text after equation[3]). As the mass resolution is increased 
the gas can cool to progressively lower temperatures. It can 


also be seen that in simulations with 10 5 and 10 6 particles 
the maximum temperature of the gas (before any AGN feed¬ 
back is initiated) converges to T ~ 10 3 ' 7 K. The reason for 
this is that at lower temperatures (T ^ 10 4 K) the gas cools 
much slower, so that there tends to be a lot of gas “piling 
up” at T ~ 10 4 K. 


3.2 Overview of numerical resolution trends 

Figure [3] shows rendered images of gas density slices through 
the simulation domain at y — 0 after 1.5 Myr. The top row 
shows the FM3h, FM4h, FM5h and FM6h runs from left 
to right respectively while the bottom row shows the FN3h, 
FN4h, FN5h and FN6h runs from left to right respectively. 
The figure clearly illustrates the increasing complexity of 
structure that can be resolved with improved resolution. In 
the low resolution runs (left panels) there is a fairly symmet¬ 
rical swept-up shell of high density gas, while the high resolu¬ 
tion runs (right panels) consist of compressed high density 
filaments and cleared low density channels through which 
hot gas can escape. 

Figure 0] shows the respective temperature slices. At 
low resolution (left panels) the feedback outflow sweeps up 
essentially everything in its path, with no cold gas left at 
small radii, whilst the hot gas is contained in the central 
regions only. In contrast, in the higher resolution runs (right 
panels) the cold gas is still present in clumps and filaments 
at small radii, whereas the heated gas escapes through low 
density channels and is now more spatially extended. Since 
it is likely that the cold gas is the source of efficient SMBH 
growth, these results show that not only SMBH feedback but 
also SMBH growth is affected by the numerical resolution 
artefacts i.e. at low resolution there is a lack of high density 
cold gas clumps. 

The simulations also show stark differences in gas ther¬ 
modynamical properties between the runs in which the feed¬ 
back is coupled to a fixed mass (the FM series of runs) versus 
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Figure 3. Density slices at y = 0 through simulation domains at 1.5 Myr. The top row shows runs FM3h, FM4h, FM5h and FM6h from 
left to right respectively while the FN3h, FN4h, FN5h and FN6h runs are shown from left to right respectively on the bottom row. 



Figure 4. Temperature slices at y = 0 through simulation domain showing gas temperature at 1.5 Myr. The top row shows runs FM3h, 
FM4h, FM5h and FM6h from left to right respectively while the FN3h, FN4h, FN5h and FN6h runs are shown from left to right 
respectively on the bottom row. 


those with a fixed neighbour number (the FN series of simu¬ 
lations). For instance, due to the differences in the feedback 
implementation, in FN6h (bottom right panels of figures [3] 
and il a factor ~ 100 times less mass is heated than in 
FM6h (top right panels of figures [3] and H} and hence the 
gas is heated to much higher temperatures in FN6h than in 
FM6h. This results in the feedback in FN6h being more ef¬ 
fective at clearing gas from the central region. However there 
is still cold, dense, in-flowing gas present at small radii. 


3.3 Impact of feedback on the ISM 

3.3.1 Resolving the ISM density structure 

Figure [5] compares the density distribution at lMyr (filled) 
and 2Myr (not filled) (i.e. before and after the AGN out¬ 
burst) for the FN6 (blue) and FN3 (red) runs. In the FN3 
runs there is very little evolution in the density distribution 
with time. However, in the FN6 runs the highest value of 
gas density reached in the simulation increases by approxi¬ 
mately two orders of magnitude, especially when IC cooling 
is included (compare the dotted and solid lines). 
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Figure 5. Comparison of gas density distributions at 1 Myr 
(filled) and 2 Myr (not filled) for FN3 (red) and FN6 (blue) runs 
with and without IC cooling (solid and dotted respectively). It is 
clear to see that in the FN3 runs the density distribution changes 
very little, while in the FN6 runs the gas can be compressed to 
considerably higher densities. 
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Figure 7. Comparison of radial velocity distributions at 1 Myr 
(filled) and 2 Myr (not filled) for FM3 (red) and FM6 (blue) 
runs with and without IC cooling (solid and dotted respectively). 
Both high and low resolution runs produce gas out-flowing with 
velocities of order 1000 km s _1 , however the high resolution runs 
(FM6) show far stronger gas inflows than the low resolution runs 
(FM3). 
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Figure 6. As in figure [5] except for the FM3 and FM6 runs. How¬ 
ever, the FM6c run exhibits a particularly high density feature 
not seen in other runs. 


Similar behaviour is seen when comparing the FM3 and 
FM6 runs in figure [6] The FM6c run exhibits a particularly 
high density feature not seen in other runs. While the mass 
in this feature corresponds to only ~ one resolution element, 
it is interesting to note its existence. The likely cause of this 
clump is two-fold. Compared with the FM6h run, the gas 
in the FM6c run can cool via IC processes allowing it to 
collapse to higher densities. Furthermore, comparing with 
the FN6 runs, the gas heated directly by the AGN feedback 
does not reach such high temperatures in the FM6 runs, 
potentially resulting in cooler clumps that can reach higher 
densities. 

In the higher resolution runs the AGN feedback is 
able to compress gas to much higher den sities, which could 
result in triggered star formation (e.g. lElbaz et~al ] hood : 



radial velocity km s 1 


Figure 8. As in figure [7] except for the FN3 and FN6 runs. Due 
to the significantly higher temperatures reached in the FN6 runs, 
compared to the FM6 runs, the gas outflows can reach much 
higher velocities. 


Gaibler et al 


Zubovas et al 


20121: iNavakshin fc Zubova3l2012l : Isilkl 1201.1 
l2013al : ICresci et, all 20151 ). The figure demon¬ 


strates that compression of the ISM into high density fea¬ 
tures is largely missed in the low resolution runs probably 
because the ISM structure is under-resolved. 


3.3.2 Resolving outflows and inflows 

The radial velocity of gas in the simulations is also affected 
by numerical resolution. Contrasting the FM3 and FM6 
runs, figure [7| shows that whilst both high and low reso¬ 
lution runs produce gas out-flowing with velocities of order 
1000 km s -1 , the same cannot be said about the in-flowing 
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gas: the high resolution runs (FM6) show far stronger gas 
inflows than the low resolution runs (FM3). The same be¬ 
haviour is found when comparing the FN3 and FN6 runs in 
figure [HI It is interesting to note that for this implementa¬ 
tion of feedback, the out-flowing gas can reach much higher 
velocities in the FN6 run than in the FN3 run (by a factor 
of ~ 10). This can be attributed to the much higher temper¬ 
atures achieved in the FN6 run. The physical reason for the 
in-flowing gas present only in the high resolution runs is the 
previously emphasised inability of the low resolution simula¬ 
tions to model the high density features properly. The high 
density clumps and filaments present at high resolution are 
artificially smoothed in lower resolution runs. This results 
in the high density gas being far less resilient to feedback in 
the low resolution runs and hence being blown away with 
the rest of the gas. Needless to say, this is a serious arte¬ 
fact as the SMBH may be fed by exactly this high density 
gas falling int o the centre of th e galaxy despite the SMBH 
feedback fe.g.. lNavakshinll2014l h 


3.4 Efficiency of feedback versus numerics 

3-4-1 The over-cooling problem 

Supernova feedback simulations show a well known 
“over-cooling problem” which af fects simulation resu lts 
and is discussed at length in I Dalia Vecchia fe Schay e 
l2012h . We, like other studies f e.g. iBooth & Schave 


20091: iM cCarth v et al.l l20ld. 2011; Ike Brun et a 1.1 l20i4 
Scha ve et al. ^OlsTTCrain et ah l201ol) . find that this prob- 


lem also exists in AGN feedback simulations. The maximum 
temperature of the gas directly heated by the feedback, that 
is the SPH neighbours of the SMBH particle in which the 
feedback energy is directly deposited, is inversely propor¬ 
tional to the total mass of the gas heated. The radiative 
cooling rate of the gas is a strong function of temperature in 
certain temperature ranges. Therefore, the impact of radia¬ 
tive cooling on the thermal evolution of this gas depends in 
a complicated fashion on the number or total mass of SPH 
particles in which the feedback energy is injected. In low 
resolution simulations it is likely that the injected energy 
is spread over an unrealistically large mass of ambient gas. 
This typically means that this feedback-heated gas cools on 
timescales much shorter than one would physically expect. 

Figure [9] shows the time evolution of the instantaneous 
maximum gas temperature for simulations in which the feed¬ 
back energy is injected into a fixed number of SPH particles 
(~ 100) during each time step (FN3 (red), FN4 (black), 
FN5 (green) and FN6 (blue)) in the top panel (a) and for 
simulations in which the feedback energy is injected into a 
fixed mass (~ 1.6 x 10 7 Mq) of gas during each time step 
(FM3 (red), FM4 (black), FM5 (green) and FM6 (blue)) 
in the bottom panel (b). The solid and dotted lines show 
runs with and without IC cooling respectively. Considering 
first the fixed number of neighbours (FN) runs in the top 
panel (a), it roughly follows that each order of magnitude 
improvement in resolution results in an order of magnitude 
increase in temperature. This is because the feedback energy 
is injected into a factor ~ 10 times less gas mass for each 
factor 10 increase in mass resolution. 

The fixed mass (FM) runs in the bottom panel (b) lead 
to more consistent results in that the maximum temperature 




Figure 9. Time evolution of the maximum gas temperature. The 
top panel (a) shows the FN3 (red), FN4 (black), FN5 (green) 
and FN6 (blue) runs whilst the bottom panel (b) shows the FM3 
(red), FM4 (black), FM5 (green) and FM6 (blue) runs. Solid and 
dotted lines indicate runs with and without IC cooling.This figure 
shows that the temperature to which gas is heated to is strongly 
dependent on the mass of gas heated. 


of gas varies much less between the different resolution sim¬ 
ulations, as is expected. However, at later times the higher 
resolution runs FM5 and FM6 do differ significantly from 
the lower resolution curves. We believe this is caused by dif¬ 
ferences in gas properties beyond the immediate feedback 
deposition region. As we know from figure 3, there is more 
dense gas near the black hole in the better resolved sim¬ 
ulations. This denser gas is hence better at radiating the 
feedback energy away than in the lower resolution runs. 

Finally we note that if included, IC cooling dominates 
the cooling function at high temperatures for gas close to 
the SMBH. This explains why the dotted and dashed curves 
in figure [9] only exhibit differences when gas is heated to high 
temperatures. 

The mode of feedback energy delivery to the ambient 
gas strongly affects the subsequent evolution of the system. 
This can be seen in figure [TO] which shows the time evolution 
of the ratio of the change in gas internal and kinetic energy 
to the energy injected by the AGN as a function of time. 
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Figure 10. Evolution of the ratio of total gas energy, E gas , to 
total AGN input energy, Fagn> from the beginning of the AGN 
outburst between 1 Myr and 2 Myr. The top panel (a) shows the 
FN3 (red), FN4 (black), FN5 (green) and FN6 (blue) runs. In this 
panel there is no apparent trend with resolution. In both panels 
the solid and dotted lines correspond to runs with and without 
IC cooling respectively. The bottom panel (b) shows the FM3 
(red), FM4 (black), FM5 (green) and FM6 (blue) runs. There is 
an evident trend that lower resolution runs retain more energy 
than high resolution runs. 


As in figure [9J in the top panel (a) the feedback energy is 
injected into a fixed number of SPH particles (~ 100) during 
each time step whereas in the bottom panel (b) the feedback 
energy is injected into a fixed mass (~ 1.6 x 10 7 M©) of gas 
during each time step. The change in total gas energy at 
time t, A E gas (t), is given as 

A E gas (t) = E gas o t ) E gas (lMyr) (9) 

where 

+ ( 10 ) 

is the sum of the kinetic and internal energy of all of the 
SPH particles at time t. The total AGN input energy, Pagn 


is given by 

pAGN = Cf^AGNGct (11) 

where Get is the time for which the AGN has been active. 
Focusing on the FN series of runs first (top panel, a), we 
find no clear trend in how much feedback energy is retained 
by the gas. The two low resolution cases, FN3 and FN4, are 
rather similar; then the higher resolution case, FN5, retains 
less energy than that, but the highest resolution case, FN6, 
retains much more energy than the low resolution cases. We 
believe that this is due to competition between two some¬ 
what oppositely directed numerical artefacts, one due to the 
over cooling problem close to the AGN and the other due to 
poor sampling of the ambient gas farther out. As the reso¬ 
lution is improved, one is able to heat a lower mass of gas 
and hence heat the gas to higher temperatures, resulting in 
longer cooling times. However at higher resolution one can 
resolve high density material which has a shorter cooling 
time than the lower density material found in the vicinity of 
the SMBH in lower resolution simulations. Given the compe¬ 
tition between the processes contributing to the gas cooling 
rate there is not necessarily a clear trend in feedback energy 
conservation with resolution. 

The fixed mass runs (FM, the bottom panel, b, of the 
figure) give a more consistent picture, with FM5 and FM6 
yielding very similar results, perhaps indicating that a de¬ 
gree of numerical convergence is taking place. Here we see 
that at higher resolution less feedback energy is retained in 
the ambient gas of the galaxy, presumably because higher 
density clumps are better resolved at higher resolutions and 
they lead to more energy being lost to radiation. These 
results show that simulations in which feedback is spread 
around a fixed mass of ambient gas should be preferred 
for numerical reasons to those where the number of AGN 
neighbours is fixed instead. However it should be noted that 
although injecting feedback energy into a fixed mass of par¬ 
ticles provides a degree of numerical convergence, it is not 
necessarily convergence towards the physically correct re¬ 
sult. A further potentially confounding factor is the energy 
imparted by the momentum of the AGN wind. This is not 
included in the models presented in this paper, however a 
purely momentum-driven wind should have a kinetic en¬ 
ergy rate E mom ~ 2a/(r]c)E wind ~ 0.01G wind . In the high- 
resolution models, this energy is comparable to the energy 
retained by the gas and may further complicate gas be¬ 
haviour. However, a more detailed investigation of the ef¬ 
fects of different AGN feedback prescriptions is beyond the 
scope of this paper. 


3-4-2 Gas ejection efficiency 

It is believed that the most important effect of AGN feed¬ 
back onto their host galaxies is to remove gas from the host 
galaxy and thus quench star formation. In this section we 
show that the ability of an AGN to clear the gas out of the 
host is greatly affected by numerical resolution at least for 
the initial conditions and parameter space investigated in 
this paper. 

To quantify the AGN feedback impact on the host, we 
first consider the evolution of the total baryonic mass en¬ 
closed within 200pc of the host’s centre. We define the frac- 
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tional change in baryonic mass as 

A7\/f< 2 oopc (^) Af< 2 oopc(^) Af< 2 oopc (lMyr) M2) 

M<2oo P c(lMyr) M< 2 oo P c(lMyr) 

where M< 2 oo P c(£) is the total baryonic mass within 200 pc, 
including gas accreted onto the black hole but not the ini¬ 
tial black hole mass ( 10 8 M©) itself. Figure fill shows the 
time evolution of AM< 2 oopc(£)/Af< 2 oo P c(lMyr) for simula¬ 
tions with 10 3 , 10 4 , 10 5 and 10 6 particles shown in black, 
blue, red and green, respectively. The fixed mass (FM, bot¬ 
tom panel, b) runs show a trend with resolution such that 
feedback becomes less effective at clearing gas out with im¬ 
proved resolution. The FM3 run which retains the largest 
fraction of energy is the most effective at clearing gas out, 
whilst the FM5 and FM6 runs, which lose the most energy, 
cannot prevent the continual infall of gas. However, if we 
consider the FN runs (top panel, a) there is no clear trend. 
Even though the FN6 run retains substantially more energy 
than the FN3 run, it is far less effective at clearing gas out. 
This can be attributed to being able to resolve structure 
in the ISM. As we improve resolution, the hot gas can es¬ 
cape through paths of least resistance leaving higher density 
clumps behind. 

Finally we attempt to further quantify the ability of 
an AGN to clear gas out by calculating the fraction of gas 
with radial velocity greater than 2cr. This is shown in figure 
da with the same colour scheme as previous plots. For the 
FM runs (bottom panel, b) there is a clear trend that as 
the resolution is improved the fraction of gas out-flowing at 
a high velocity becomes smaller. However for the FN runs 
(top panel, a), this trend breaks down for the highest resolu¬ 
tion run, FN6. This suggests that the cooling as well as the 
ability to resolve structure plays an important role in deter¬ 
mining how effectively gas can be blown out of the galaxy. 
As discussed when considering figures [H and [9] with respect 
to the FN5 and FN6 runs, the gas in the FN6 runs is heated 
to higher temperatures than that in the FN5 runs. The cool¬ 
ing in the FN6 runs is less efficient and thus the hot feedback 
bubble conserves more energy than in the FN5 runs and is 
able to drive more powerful outflows. 


4 DISCUSSION AND SUMMARY 

We have studied the effect of AGN feedback on a multi¬ 
phase interstellar medium and how such feedback is affected 
by numerical resolution. The resolution has two competing 
effects on the results. Unsurprisingly, the density structure 
is better resolved at higher resolutions, so that there is more 
hot low-density and also cold high-density gas than in low 
resolution simulations. Low resolution simulations thus tend 
to expel all the gas from the centres of host galaxies, whereas 
higher resolution simulations do not; they also show dense 
clumps that are more resilient to feedback and may continue 
to fall inward while most of the gas is driven out. On the 
other hand, the over-cooling problem may affect the gas in 
the immediate vicinity of the SMBH, actually reducing feed¬ 
back efficiency in low resolution simulations. We now discuss 
these opposing effects in detail. 



time (Myrs) 



Figure 11. Time evolution of the fractional mass change within 
the central 200 pc of the system. The top panel (a) shows the 
FN3 (red), FN4 (black), FN5 (green) and FN6 (blue) runs. As in 
figure [Tol there is no apparent trend with resolution, however, like 
in the FM runs, the FN3 and FN4 runs are the most efficient at 
removing gas. The bottom panel (b) shows the FM3 (red), FM4 
(black), FM5 (green) and FM6 (blue) runs. As the resolution 
is degraded the feedback becomes more efficient at clearing gas 
from the central regions. In both panels the solid and dotted lines 
correspond to runs with and without IC cooling respectively. 


4.1 Resolving the multiphase ISM and outflow 
properties 

The inferred properties of outflows in our simulations 
strongly depend upon the resolution and how the feedback 
energy is coupled to the ISM. At low resolution the feed¬ 
back sweeps up all the material in SMBH vicinity into an 
outflow with a modest velocity (~ 1000 km s -1 ). In con¬ 
trast, at high resolutions only some of the neighbouring gas 
is launched into the outflow. The outflows can reach higher 
radial velocities (up to ~ 5000 km s -1 ); however, cold dense 
filaments may continue to fall in and feed the SMBH. Higher 
resolution SPH simulations thus lead to a much more var¬ 
ied outcome for the gas in the galaxy in every measurable 
quantity, e.g., density, temperature and velocity. This is po¬ 
tentially very important for SMBH growth since even a tiny 
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Figure 12. Time evolution of the fraction of gas with radial ve¬ 
locity greater than 2 cr. The top panel (a) shows the FN3 (red), 
FN4 (black), FN5 (green) and FN6 (blue) runs. Again, as in pre¬ 
vious figures there is no apparent trend with resolution with both 
the FN3 and FN6 runs producing outflows containing similar frac¬ 
tions of gas. The bottom panel (b) shows the FM3 (red), FM4 
(black), FM5 (green) and FM6 (blue) runs. As in figure fill there 
is an apparent trend with resolution, with a greater fraction of 
gas being out-flowing at a high velocity in the lower resolution 
runs than in the high resolution runs. In both panels the solid 
and dotted lines correspond to runs with and without IC cooling 
respectively. 


mass of ambient gas is sufficient to increase the SMBH mass 
significantly. 

In the same vein, cosmological simulations often have 
to invoke high stellar feedback e fficiencies in order to pro¬ 
duce observed galactic winds ([Schave et al.l 120151 ). One 
possible reason for this, in addition to the over-cooling 
problem, is that low resolution inhibits low density chan¬ 
nels through which outflows can escape to reach galac¬ 
tic scales. Alternatively, in some cosmological and galaxy- 
scale simulations, such winds are hydrodynamically decou¬ 
pled from t he ISM so that they can freely stream to galac- 
tic scales (Springel & Hernquist 2003; Oppenheimer & Dave 

120061 : lOppenheimer et alJl20ld : IPuchwein fe Springelll2013f h 

While such ad-hoc prescriptions allow one to produce real¬ 


istic outflows at large radii, one loses any information re¬ 
garding the direct interaction of feedback with the ISM. It 
is therefore clear that cosmological simulations are unable 
to provide detailed insights into the feedback mechanisms 
themselves and the best that one can hope for is that their 
effects on resolvable scales are modelled correctly. 

4.2 Cooling of the feedback bubble 

It is evident that the ability of the feedback bubble to con¬ 
serve its energy has a significant impact on how efficient it 
is in destroying the host galaxy. This can depend intimately 
both upon the temperature to which gas is heated directly 
by the AGN and the processes through which the feedback 
heated bubble cools. Considering the first point, we have 
shown that there are stark differences in the properties of 
the feedback depending upon the mass of the gas heated by 
AGN feedback. This leads us to the question, what tempera¬ 
ture is correct? Analytical theory suggests that wind shocks 
have temperatures of 10 10 — 10 11 K, however, as illustrated 
by figure [9] in order to reach such high gas temperatures in 
a cosmological simulation, a gas mass of £ 1 SPH particle 
would need to be heated. Therefore, in such simulations, gas 
is typically heated to lower temperatures, which could po¬ 
tentially result in incorrect energy conservation in the hot 
bubble. Such problems may be mitigated to some degree by 
tuning the efficiency of the feedback or artificially turning 
off radiative cooling in order to match observations. Whilst 
such procedures can lead to correct large scale properties of 
the galaxies, e.g. correct stellar masses and the Mbh — cr 
relation, one clearly loses predictive power if observations 
need to be used to calibrate the models. 

With regard to our second point, i.e. the relevance of 
different cooling mechanisms, the inclusion or absence of IC 
processes can have a big effect. Comparing the FN 6 c and 
FN 6 h runs, it is clear from figures flQl fill and fl2l that the 
hot bubble retains more energy and clears out more gas 
when IC cooling is neglected. It is therefore important to 
understand which scenario is mor e physically motivated. 
iFaucher-Giguere & Quataertl (|2Q12h have shown that given 
the high temperature and low density properties of shocked 
outflows, the electrons and ions are thermally decoupled. 
The electron cooling timescale is shorter than the electron- 
ion thermalisation timescale and therefore it is the latter 
that determines the cooling rate of the gas, with IC cooling 
becoming ineffective. This suggests that the runs in which 
IC cooling is neglected are more physically motivated. How¬ 
ever, even if the mass and temperature of the feedback bub¬ 
ble matches those expected from analytical theory, the sim¬ 
ulated feedback bubble likely has a higher density than the 
actual shocked wind bubble would, because we are directly 
heating ISM gas. Therefore the cooling of the hot bubble 
and its interaction with the ambient ISM may still be in¬ 
correctly modelled. Thus we conclude that direct physically 
self-consistent modelling of AGN feedback heating and cool¬ 
ing on small scales is still beyond the reach of modern nu¬ 
merical capabilities. 


4.3 Star formation during an AGN outburst 

Figures [5] and [ 6 ] clearly show that gas can be compressed to 
high densities by an AGN outflow. The presence of dense gas 
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can result in additional star formation, which both quantita¬ 
tively and qualitatively changes the properties of the AGN 
host galaxy. Similar aspects of AGN-triggered star formation 
hav e already been explored bv lNavakshin fc Zubovasl (|2012h 
and IZubovas et al.l (l2013aT ) , who found that significant star 
formation can occur both in the cooling out-flowing medium 
and in the compressed disc of the host galaxy. Our results 
show that any density contrasts can be enhanced by AGN 
outflows. 

It is important that one needs simulations with suffi¬ 
cient resolution to recover the compression effect in numer¬ 
ical simulations (Figures [5] and [6j) . Large-scale simulations 
with low numerical resolution typically miss this effect and 
hence over-predict the negative (gas removal) effect of AGN 
feedback. Even in high-resolution simulations, star forma¬ 
tion in dense clumps is difficult to track during the AGN 
outflow if one employs a heatin g-cooling prescription such 
as that of ISazonov et al.l (120051 ). which includes Compton 
heating. This prescription assumes that gas is optically thin, 
which is a good approximation for the low density ISM, but 
not for the dense clumps. As a result, the clump tempera¬ 
ture, in general, stays too high (i.e. above the temperature 
set by equation [3]) and fragmentation is slower than it would 
be with a proper radiative transfer treatment. 

The lack of AGN-triggered star formation in low- 
resolution simulations of galaxy evolution presents two chal¬ 
lenges: quantitative and qualitative. The quantitative chal¬ 
lenge is the issue of reconciling the star formation histo¬ 
ries of simulated and observed galaxies. If one channel of 
triggered star formation is missed in simulations, the other 
star formation channels have to be proportionately enhanced 
(for example, by adopting higher star formation efficien¬ 
cies or lower density thresholds) in order to reproduce the 
galaxy stellar mass functions of present-day galaxies. The 
qualitative challenge is arguably more important: AGN out¬ 
flows create dense star-forming gas where there was none, 
i.e. affect the location of star-forming regions in the galaxy. 
This process directly affects the morp hology of the starburst 
and the dynamics of new-b orn stars dNavakshin fe Zubovasl 
120121 : IZubovas et al.ll2013bh . Both of these effects are missed 
in low-resolution simulations; however, they can be used as 
strong indicators of positive AGN feedback. 

One region where AGN-induced star formation may 
be particularly important is galaxy centres. T hese regions 
typically contain dense gas discs or rings (e.g. iBoker et al.l 
12008I L which often show clumpy structures and embedded 
young star clusters. It is generally accepted that star forma¬ 
tion in these regions is induced by shocks caused by mat¬ 
ter in-falling via galactic bars from larger radii. However, 
the presence of young clusters and the lack of azimuthal 
age g radient in some o f these systems complicate this pic¬ 
ture ([Boker et al.ll2008h . Another trigger of star formation 
in these systems could be AGN outflows. As these outflows 
expand perpendicularly to the dis c plane due to lower den¬ 
sity ([Zubovas &; Navakshinl l2014h . they significantly com¬ 
press the gas in the midplane; in addition, ram pressure 
of the AGN wind pushes the disc gas into a narrower ring, 
which is more prone to gravitational instability (Zubovas, 
submitted). In this way, the density contrast between the 
disc and its surroundings is also enhanced, much like the 
density contrast between different regions in the simulations 
presented in this paper. 


4.4 Black hole growth and the Mbh — cr relation 

The general trend in the simulations we present is that at 
higher resolution less gas is cleared out than at lower reso¬ 
lutions. In low resolution runs the outflow sweeps up every¬ 
thing in its path, creating a sharp cut-off radius between 
out-flowing material and in-flowing material. However in 
high resolution runs the outflow only sweeps up low den¬ 
sity material, whilst high density material can continue to 
flow inward. This could lead to very different feeding cycles 
for the black hole. The supply of material to the black hole 
is completely cut off and cleared to large radii in the low 
resolution runs, whereas in the high resolution runs clumps 
and filaments can remain in-flowing at small radii and thus 
continuously feed the black hole. This sets up a scenario 
in which feedback is “all or nothing” at low resolution but 
more diluted at high resolution, with feeding becoming in¬ 
terminable up to the point that the gas can form stars. 

In the high resolution scenario, the high density clumps 
will only be acted upon b y the momentum of the AGN wind 
(BNH14. lNavakshir3l2014h with the energy escaping through 
low density channels. By requiring the ram pressure of the 
AGN outflow to exceed the gravitational force of the bulge 
acting on all of the clumps along the line of sight from a 
SMBH and setting a maximum threshold density for the 
clumps (assum ed to be the dens ity at which they undergo 
star formation) iNavakshinl (I2Q14T ) finds a critical black hole 
mass in order to clear out the cold gas of 

Merit ~ 2.2 x 10 s M 0 4 oo (13) 

comparable to the o bserved Mbh — cr relations (e.g., 
iKormendv fe Holl2013h . 


4.5 Comparison with previous work 

Recent work (| Wagner et al.l [20131 . BNH14) has shown that 
the structure of the ISM can impact upon the ability of 
an AGN to clear out gas and hence quench star formation. 
BNH14 present high-resolution simulations of an UFO im¬ 
pacting upon an inhomogeneous, turbulent medium and find 
new processes such as energy leakage and separation of en¬ 
ergy and mass flows within the ISM. The shocked outflows 
escape via paths of least resistance, leaving the high density 
gas, which is difficult to expel, largely intact. Such processes 
have previously been missed in analytical models and cosmo¬ 
logical simulations mainly because the multiphase nature of 
the ISM is not resolved and has to b e implemented as a sub - 
resolution model. For example ISoringel fe Hernauistl (12003 4 
include a sub -gr id mu ltiphase model for star formation while 
iMurante et al.l (120101 ) include a non-equilibrium model that 
includes the three ISM phases for each SPH particle. Un¬ 
fortunately such methods mean that the intricate structure 
one expects is washed out due to low resolution. 

This work builds on that of BNH14 by implement¬ 
ing a continuous, Eddington limited feedback outburst, 
rather than a single hot bubble. Furthermore, we have 
explored the role of IC cooling against the AGN ra¬ 
diation field, which has been highlighted in the liter¬ 
ature as a key feature in understanding the impact 
of UFOs ([Faucher-Giguere fe Quataertl l2012h. This work 
adds to th e growing body o f work (e.g. I Wagner et ahl 
120121 . 120131 : ICosta et all 1 2014 iGabor & Bournaudl 20 14 
























































The resolution bias 13 


IZnhovas fc Na,vakshinlEoi 4. BNH14) highlighting that the 
AGN environment can be just as important as the AGN 
feedback mechanism itself when modelling galaxy evolution. 
In BNH14, the main aim was to understand the physics of 
the interaction of an outflow with the multiphase ISM, how¬ 
ever in this work we have focused more on how resolution 
can affect this interaction. 


The role of the ISM and its impact on AGN feedback 
has been studied by a n umber of au t hors b oth for feedback 
in the form of jets ( e.g. lWagner e t al.l2012h a nd UF Os (e.g. 
I Wagner et alJl20l3 . BNH14j. IWagner et al.l (l2012l l present 
high resolution simulations of jet feedback in a clumpy ISM. 
They found that if the volume filling factor of the clouds is 
less than 0.1 then the hot feedback bubble can expand as in 
the energy driven limit. Clouds smaller than ~ 25 pc are de¬ 
stroyed and dispersed, leading them to argue that feedback 
prescriptions in cosmological simulations should provide a 
good description of this regime as a source of negative feed¬ 
back. However if clouds are larger than ~ 25 pc they are 
more resilient to the feedback. In agreement with this work 
they find that the clouds can be compressed, potentially 
triggering star formation. Suc h behaviou r is misse d in cos¬ 
mological simulations. Whilst I Wagner et al.l (2012) suggest 
a physical setup in which feedback prescriptions in cosmo¬ 
logical simulations may produce correct results, we provide 
a direct comparison of the nature of feedback when simu¬ 
lated at low resolution (similar to cosmological simulations) 
and up to three orders of magnitude higher resolution. We 
have found that across such a resolution range there are 
marked differences in the evolution of the feedback, caused 
by a combination of effects including the ability to resolve 
structure and the thermal and physical properties of the hot 
feedback bubble. 


Further work on scales simulating whole galaxies 
has shown that large s cale structure such as a disk 
([Gabor fe Bournaudll2014h or filaments ([Gosta et al.ll2014h 
can also reduce the ability of AGN feedback to remove gas 
from the galaxy. Such structure should be resolved in cosmo¬ 
logical simulations, however resolution effects can still im¬ 
pact upon the properties of the feedback. As we have shown, 
the temperature to which gas is directly heated by feedback 
as well as its density can affect the cooling and thus the effi¬ 
ciency of the feedback. Given such large differences between 
the physical properties of hot feedback bubbles in cosmolog¬ 
ical simulations and those expected in reality we should pose 
the question when, if ever, such processes can be included 
in these simulations. Following a Moore’s law approach the 
number of particles used in cosmological N-body simulations 
approximately doubles every 16 months. This would suggest 
that an increase in the mass resolution by 3 orders of mag¬ 
nitude could be achieved in ~ 13 years. However, given the 
fact that algorithms typically scale worse than O(N), and 
also considering that the silicon chip capacity is limited, this 
is an extremely optimistic estimate. It is therefore likely that 
such an improvement in resolution would take much longer 
to achieve and depends upon the efficiency with which sim¬ 
ulators and programmers can harness the power of parallel 
processing and other technological advances. 


4.6 Implications for cosmological simulations 

A caveat to the results presented in this paper is that our 
simulations do not include self-regulation of the AGN feed¬ 
back, which plays an important role in galaxy evolution. 
Instead our simulations only model a single, 1 Myr long, Ed¬ 
dington limited AGN feedback event which is not linked to 
the gas content of the host galaxy. It may therefore be argued 
that our results on the numerical artefacts in AGN feedback 
efficiency do not have direct implications for cosmological 
simulations. In such simulations, the system will undergo 
multiple feedback events over cosmological timescales. The 
rate at which a black hole injects energy into the host galaxy 
ISM is coupled to the gas accretion rate rh ac cr through the 
equation 


E = £bh(Mbh ,rrt accr)hlaccr C 


(14) 


where cbh(Mbh, rh acC r) is as defined in the introduction 
and can be a funct ion of Mbh and rh ac cr- For example 
ID avis &; Laorl (l201lh determined e r in 80 quasars by using 
their bolometric luminosities and absolute accretion rates, 
calculated using thin accretion disk model spectral fits, find¬ 
ing a scaling with Mbh such that e r = 0.089Ms' 52 , where Mg 
is Mbh in units of 10 8 Mq. Often, however, cbh(Mbh, ffiaccr) 
is set to be a constant for simplicity. 

The coupling of the AGN feedback to the gas content 
of the galaxy through equation [14] leads to self-regulation of 
the SMBH growth and feedback resulting in the correct E 
such that the feedback driven outflows balance mass inflow. 
This, therefore, does not uniquely establish Mbh but rather 
the product cbhMbh (because m acc r is usually limited to 
the Eddington accretion rate). To reproduce the observed 
black hole correlation s , one fixes the value of cbh (e.g., 
iBooth fe Schavel 120091 . l20ld : ISchave et al.l l2015l j. Further, 
provided that cbh is set to a value within a suitable range, 
the observed SMBH scaling relations can be reproduced 
despite large d ifferences in resolution and sub-grid pre¬ 
scriptions (e.g., iDi Matt eo et al .l I2QQ5I: ISp ringel et al J |2QQ5|: 


Siiac ki et all 2007 ; Booth &; Schavel2009f Schave et al.l 2Q1C , 


2015), although some fine tuning may be required (see dis¬ 
cussion in the introduction). A key element of self-regulation 
is that the phys ical properties of the galaxies, such as the 
stellar mass (e.g. jDi Matteo et alJl2005l:(Springel et alJl2005l : 

I Sii acki et al.l 120071 : Booth fe Schavel 20091 . 201 Oh or AGN 
driven outflow rates (ISchave et alJl2015h do not depend upon 
the chosen value of cbh- The result of this is that any de¬ 
pendencies that cbh has on resolution whould not effect the 
global properties of the galaxy due to self-regulation, al¬ 
though may lead to changes in the SMBH mass. 

As an alternative to tuning efficiencies, a number of 
authors have attempted a more physically constrained ap¬ 
proach to AGN feedback, in which cbh(Mbh, rh acC r) is not a 
fre e par a meter . For example, we here followed the model 
of I Kind (l2005l l in setting €f — 0.05. Additionally, there 
is growing evidence that AGN should u ndergo separate 


quasar and radio modes of feedback (e.g., IChurazov et al 


2005; iHeinz et al. 1 120051 : ICroton et al.l [2006 ; Ishi bashi et al 


20141 5 depending on the Eddington ratio, ra acC r/MEdd, 


each with differing values for Cf. It has also been sug¬ 


gested that e r depends upon rh accT (e.g.. iNaravan fe Yi 
1995:lMahadevanlll997l:ICiotti fe Ostriker 2001; Ci otti et al 


20091 5, with some cosmological simulations already attempt- 
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ing t o include this additional physic s (e.g., ISiiacki et al.l 


pnysic s 

l l2014h . 


Further, an effect 


2007 . 120151 : IVogelsberger et ahl 
that is not typically taken into account in galaxy forma¬ 
tion simulations is that of the black hole spin, which can 
lead to variations in the radiative efficiency in the range of 
0.055 < e r < 0.42. 

We believe that the future of the field is in these 
more physically motivated approaches which would hope¬ 
fully provide predictions constraining the physics of black 
hole growth and feedback. Despite our simulations lack¬ 
ing self-regulated AGN feedback our results are still im¬ 
portant for such an approach. As discussed in the in¬ 
troduction some fine tuning of sub-grid feedback parame¬ 
ters can still be necessary when cha nging resolution (e.g., 
ISchave e t al. 20151: ICrai n et al.l [201 5l ). sub-grid fSM models 
(e.g.. iBooth & - Schaye 200(1) or cooling prescriptions (e.g., 
ISiiacki et al. 20151 ). When comparing results of simulations 
with different AGN feedback physics to the observations, 
one must be acutely aware of numerical artefacts that may 
skew the interpretation of such comparisons. 

Further we note the potential dependence that cbh has 
with the spatial resolution of the SMBH surroundings, tf a 
simulation probes this parameter on sub-pc scales, then the 
factor will determine the efficiency of BH wind production; 
on pc scales, the factor tells us something about the cou¬ 
pling between the wind and the surrounding fSM (perhaps 
about the dumpiness of the fSM); on scales of tens or hun¬ 
dreds of pc, the factor also encompasses the thermal effects 
(mostly cooling, but perhaps also heating by the AGN ra¬ 
diation field) of the gas surrounding the AGN. Therefore 
simulations with different spatial resolution might be prob¬ 
ing different processes which contribute to AGN feedback 
efficiency. This is an important point to consider when in¬ 
terpreting constrained values of cbh. 


5 SUMMARY 

tn this paper we have studied the effect of an Eddington- 
limited AGN outburst on a multiphase turbulent fSM, with 
particular focus on the effects of numerical resolution, tn 
general, at higher numerical resolution, more dense clumps 
and also voids through which the feedback can escape are 
found. This reduces the efficiency with which AGN feedback 
clears out the host’s gas. At low resolution this behaviour 
is lost as the feedback sweeps up essentially all the gas in 
its path. Additionally, depending on uncertain physical de¬ 
tail of the radiative cooling function for the gas heated by 
AGN feedback, numerical resolution also affects the amount 
of AGN feedback energy lost to radiation, and it is not pos¬ 
sible to say whether it will increase or decrease the feedback 
efficiency in a general case. It is therefore plausible that reso¬ 
lution dependent effects alter the efficiency of AGN feedback 
in such a way that it is difficult to attach solid physical mean¬ 
ing to constrained values of cbh = CfC r . We also note that 
although over cosmological timescales self-regulation results 
in consistent galaxy properties and outflow rates irrespective 
of the chosen feedback efficiency, our simulations illustrate 
certain physical processes, such as energy leakage through a 
clumpy ISM, that can only be modell ed at sufficiently hig h 
resolution. Finally, in agreement with ISchave et ahl ((20151 ). 
we therefore suggest caution when trying to “invert” the 


results of cosmological simulations (usually tuned to fit ob¬ 
servations) to learn about certain physical aspects of AGN 
feeding and feedback. 
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